;+
; NAME:
;       MAKE_SOURCE_HISTS
;
; PURPOSE:
;       Creates the source histograms needed for the EMAF paper.
;
; CATEGORY:
;       EMAF Plotting.
;
; CALLING SEQUENCE:
;       MAKE_SOURCE_HISTS
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       FCOLOR -- Name of the color to be used for FAR KDIST objects
;                 in source_hists3.eps.
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE  [Generates .eps files]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  05/01/12, TPEB -- Initial version.
;       Modified: 11/14/12, TPEB -- Moved to emaf/ directory and make
;                                   compliant.
;       Modified: 01/10/13. TPEB -- Modified Near/Far PostScript to
;                                   mirror the panels in the overall
;                                   one, adding in GLAT (as per
;                                   request of J. Bally).  Cleaned up
;                                   code so I can figure out what is
;                                   done where!
;       Modified: 03/12/13, TPEB -- Moved to new OMNI_ code.
;
;-

PRO MAKE_SOURCE_HISTS, FCOLOR=fcolor
  
  IF n_elements(fcolor) EQ 0 THEN fcolor='BLK5'
  ncolor = 'Opposite'
  IF n_elements(s) EQ 0 THEN $
     s = read_bgps_csv('./obsolete/bgps_distance_database.csv',/ver)
  IF n_elements(emaf) EQ 0 THEN restore,'./local/BGPS_V1_emaf.sav',/ver
  n = n_elements(emaf)
  restore,'./local/BGPS_V1_pvec.sav',/ver
  eind = where(constrain.emaf)
  pvec = pvec[eind]
  constrain = constrain[eind]
  
  ;; Limit Bolocat to the GLIMPSE overlap region for comparison
  ind = WHERE(s.glon_peak GE 7.5 AND s.glat_peak LE 65.25 AND $
              abs(s.glat_peak) LE 1.0, ngl)
  print,'Ngl: ',ngl
  s = s[ind]
  
  
  ;;================
  ;; General plotting defaults
  cs = 0.7
  xgap = 0.034
  bls = 2
  mls = 3
  linsize=0.5
  panel = 0.9
  
  ;;==========================================================================
  ;; Plot the Total Sample histograms
  
  
  myps,'./emaf_paper/plots/source_hists_f6a.eps',xsize=2.1,ysize=2.1,/bw
  ;; multiplot_xm,[3,1],xgap=xgap,mpcharsize=cs*0.8
  
  
  ;;====================================
  ;; Panel 1: GLAT Distribution
  bin = 0.1
  sbin = string(bin,format="(F0.1)")
  
  ;; Plot sample GLAT distribution
  plothist,emaf.b,bin=bin,xtit='Galactic Latitude [deg]',yr=[0,140],/yst,$
           ytit='N per '+sbin+' deg bin',charsize=cs,thick=3,xr=[-1.05,1.05],$
           xmargin=[7,1],ymargin=[4,1]
  
  ;; Get full Bolocat distro, plot at 1/10 value
  plothist,s.glat_peak,xarr2,yarr2,/noplot,bin=bin
  
  hdx = mean(xarr2[1:*]-xarr2[0,*])/2.
  
  FOR kk=0,n_elements(xarr2)-1 DO BEGIN
     xpol = [xarr2[kk]-hdx,xarr2[kk]-hdx,xarr2[kk]+hdx,xarr2[kk]+hdx] > $
            !x.crange[0] < !x.crange[1]
     ypol = [0,yarr2[kk],yarr2[kk],0]/10.
     cgPolygon,xpol,ypol,fcolor='BLK5',/fill,color='BLK5',thick=3
  ENDFOR
  
  ;; Replot Sample, to get on top
  plothist,emaf.b,bin=bin,/over,thick=3
  
  ;; Replot bottom x-axis to be seen atop the filled histogram
  cgAxis,xaxis=0,xtickformat='blank_axis',/xst
  ;; cgAxis,yaxis=0,ytickformat='blank_axis',/yst
  ;; cgAxis,yaxis=1,ytickformat='blank_axis',/yst
  
  
  ;; Legend and panel label
  ;; al_legend,/top,/right,box=0,['Bolocat x 0.1'],$
  ;;           charsize=cs*0.85,thick=3,color='BLK5',psym=15,symsize=2.0
  al_legend,/top,/right,box=0,['(a)'],charsize=cs,/clear
  
  
  myps,/done
  myps,'./emaf_paper/plots/source_hists_f6b.eps',xsize=2.1,ysize=2.1,/bw

  
  
  ;;====================================
  ;; Panel 2: BGPS Flux
  bin = 0.1
  sbin = string(bin,format="(F0.1)")
  
  plothist,alog10(emaf.s_ap),xarr,yarr,/noplot,bin=bin,$
           color=cgColor('Opposite')
  logr = 10.^[-1.45,max(xarr)+0.1]
  
  cgPlot,10.^(xarr),yarr,/NODATA,xtit='BGPS 40" Flux Density [Jy]',/xlog,$
         ytit='N per log(S)='+sbin+' bin',xr=logr,/xst,yr=[0,140],/yst,$
         charsize=cs,xmargin=[7,1],ymargin=[4,1]
  
  ;; Make filled histogram for the FAR sources
  plothist,alog10(s.flux_40),xfarr,yfarr,/noplot,bin=bin
  nbin = n_elements(xfarr)
  dx = mean(xfarr[1:*]-xfarr[0:*])
  hdx = dx /2.d
  
  ;; Plot this the hard way, box by box:
  FOR kk=0,nbin-1 DO BEGIN
     x = 10.^( [xfarr[kk]-hdx,xfarr[kk]-hdx,xfarr[kk]+hdx,xfarr[kk]+hdx] < $
               !x.crange[1] > !x.crange[0] )
     y = [0,yfarr[kk],yfarr[kk],0]/10.
     cgPolygon,x,y,/fill,color=fcolor,thick=3
  ENDFOR
  
  ;; Replot bottom x-axis to be seen atop the filled histogram
  cgAxis,xaxis=0,xtickformat='blank_axis',/xst,/xlog
  cgAxis,yaxis=0,ytickformat='blank_axis',/yst
  cgAxis,yaxis=1,ytickformat='blank_axis',/yst
  
  ;; Replot EMAF histogram atop...
  yarr = [0,yarr,0]
  dx = mean(xarr[1:*]-xarr[0:*])
  xarr = [min(xarr)-dx,xarr,max(xarr)+dx]
  cgOplot,10.^(xarr),yarr,psym=10,thick=3
  
  ;; Vertical line
  ave = mean(emaf.s_ap)
  med = median(emaf.s_ap)
  vline,med,linestyle=mls,thick=3
  
  
  ;; Legend and panel label
  ;; al_legend,/top,/right,box=0,['Median: '+string(med,format="(F0.3)")+' Jy',$
  ;;                              'Bolocat x 0.1'],$
  ;;           charsize=cs*0.85,linestyle=[mls,bls],thick=3,psym=[0,15],$
  ;;           symsize=2.,color=['Black','BLK5'],linsize=linsize
  al_legend,/top,/right,box=0,['(b)'],charsize=cs,/clear
  
  
  myps,/done
  myps,'./emaf_paper/plots/source_hists_f6c.eps',xsize=2.1,ysize=2.1,/bw
  
  
  ;;====================================
  ;; Panel 3: Mid-Infrared Contrast
  bin = 0.02
  sbin = string(bin,format="(F0.2)")
  
  plothist,emaf.c,bin=bin,color='Black',xr=[0,0.75],charsize=cs,thick=3,$
           xtit='Mid-Infrared Contrast',ytit='N per '+sbin+' contrast bin',$
           yr=[0,85],/yst,xmargin=[7,1],ymargin=[4,1]
  
  ave = mean(emaf.c)
  med = median(emaf.c)
  vline,med,linestyle=mls,thick=3
  
  ;; al_legend,/top,/right,box=0,['Median: '+string(med,format="(F0.3)")],$
  ;;           charsize=cs*0.85,linestyle=[mls],thick=3,linsize=linsize
  al_legend,/top,/right,box=0,['(c)'],charsize=cs,/clear
  
  print,'CONTRAST: ',m4_stat(emaf.c)
  
  myps,/done
  
  
  
  
  
  
  
  
  ;;====================================================
  ;; Near / Far plots
  
  ;; Restrict consideration to "well-constrained" sources...
  cind = where(constrain.post)
  pvec = pvec[cind]
  constrain = constrain[cind]
  emaf = emaf[cind]
  
  near = pvec.stat.kdar EQ 'N'
  far  = pvec.stat.kdar EQ 'F'
  tan  = pvec.stat.kdar EQ 'T'
  
  
  nind = where(near,nn)
  find = where(far+tan,nf)
  tind = where(tan,nt)
  print,'NEAR: ',nn,'  FAR: ',nf,'  TAN: ',nt
  

  myps,'./emaf_paper/plots/source_hists3_f8a.eps',xsize=2.1,ysize=2.1,/bw
  ;; myps,'./emaf_paper/plots/source_hists3.eps',xsize=15,/bw
  ;; multiplot,[3,1],xgap=xgap
  
  
  ;;====================================
  ;; Panel 4: GLAT -- NEAR vs FAR
  bin = 0.1
  sbin = string(bin,format="(F0.1)")
  
  plothist,emaf[nind].b,bin=bin,color=ncolor,xr=[-1.05,1.05],$
           xtit='Galactic Latitude [deg]',ytit='N per '+sbin+' deg bin',$
           charsize=cs,thick=3,yr=[0,90],/yst,xmargin=[7,1],ymargin=[4,1]
  
  ;; Far sources, filled histogram w/ median line
  plothist,emaf[find].b,fx,fy,bin=bin,color=fcolor,thick=3,/overplot,/fill,$
           fcolor=fcolor
  vline,median(emaf[find].b),linestyle=4,thick=3,color=fcolor
  
  ;; fyf = mpfitpeak(fx,fy,Af,nterms=3)
  ;; cgOplot,fx,fyf,thick=3,color='cyan'
  
  ;; Near sources, outline histogram w/ median line
  plothist,emaf[nind].b,nx,ny,bin=bin,color=ncolor,thick=3,/overplot
  vline,median(emaf[nind].b),linestyle=3,thick=3,color=ncolor
  
  ;; nyf = mpfitpeak(nx,ny,An,nterms=3)
  ;; cgOplot,nx,nyf,thick=3,color='red'
  ;; al_legend,/bottom,/right,box=0,textcolor=['red','cyan'],$
  ;;           [string(An[1],An[2],format="(F0.2,' ',F0.2)"),$
  ;;            string(Af[1],Af[2],format="(F0.2,' ',F0.2)")]
  
  ;; Replot bottom x-axis to be seen atop the filled histogram
  cgAxis,xaxis=0,xtickformat='blank_axis',/xst
  ;; cgAxis,yaxis=0,ytickformat='blank_axis',/yst
  ;; cgAxis,yaxis=1,ytickformat='blank_axis',/yst
  
  print,'CONTRAST: ',median(emaf[nind].c),median(emaf[find].c)
  
  ;; Legend and panel label
  al_legend,/top,/right,box=0,linestyle=[0,0],thick=4,color=[ncolor,fcolor],$
            ['Near','Far + Tan'],charsize=cs*0.85,linsize=linsize
  al_legend,/top,/left,box=0,['(a)'],charsize=cs,/clear
  
  myps,/done
  myps,'./emaf_paper/plots/source_hists3_f8b.eps',xsize=2.1,ysize=2.1,/bw
  


  ;;====================================
  ;; Panel 5: BGPS Flux -- NEAR vs FAR
  bin = 0.1
  sbin = string(bin,format="(F0.1)")
  
  plothist,alog10(emaf[nind].s_ap),xarr,yarr,/noplot,bin=bin,$
           color=cgColor('Opposite')
  yarr = [0,yarr,0]
  dx = mean(xarr[1:*]-xarr[0:*])
  xarr = [min(xarr)-dx,xarr,max(xarr)+dx]
  
  logr = 10.^[-1.45,max(xarr)+0.1]
  
  ;; Set up the plot environment
  cgPlot,10.^(xarr),yarr,/NODATA,xtit='BGPS 40" Flux Density [Jy]',/xlog,$
         ytit='N per log(S)='+sbin+' bin',xr=logr,/xst,yr=[0,100],/yst,$
         charsize=cs,xmargin=[7,1],ymargin=[4,1]
  
  ;; Make filled histogram for the FAR sources
  plothist,alog10(emaf[find].s_ap),xfarr,yfarr,/noplot,bin=bin
  nbin = n_elements(xfarr)
  dx = mean(xfarr[1:*]-xfarr[0:*])
  hdx = dx /2.d
  
  ;; Plot this the hard way, box by box:
  FOR kk=0,nbin-1 DO BEGIN
     
     x = 10.^( [xfarr[kk]-hdx,xfarr[kk]-hdx,xfarr[kk]+hdx,xfarr[kk]+hdx] < $
               !x.crange[1] )
     y = [0,yfarr[kk],yfarr[kk],0]
     cgPolygon,x,y,/fill,color=fcolor,thick=3
  ENDFOR
  vline,median(emaf[find].s_ap),linestyle=4,thick=3,color=fcolor
  
  
  ;; Make open histogram for NEAR sources
  cgOplot,10.^(xarr),yarr,psym=10,thick=3,color=ncolor
  vline,median(emaf[nind].s_ap),linestyle=3,thick=3,color=ncolor
  
  print,'FLUX DENSITY: ',median(emaf[nind].s_ap),$
        median(emaf[find].s_ap)
  
  ;; Replot bottom x-axis to be seen atop the filled histogram
  cgAxis,xaxis=0,xr=logr,/xst,xtickformat='blank_axis',xthick=3
  cgAxis,yaxis=1,xr=logr,/yst,ytickformat='blank_axis',ythick=3
  
  ;; Legend and panel label
  al_legend,/top,/right,box=0,linestyle=[0,0],thick=4,color=[ncolor,fcolor],$
            ['Near','Far + Tan'],charsize=cs*0.85,linsize=linsize
  al_legend,/top,/left,box=0,['(b)'],charsize=cs,/clear

  myps,/done
  myps,'./emaf_paper/plots/source_hists3_f8c.eps',xsize=2.1,ysize=2.1,/bw
  
  
  ;;====================================
  ;; Panel 6: Mid-Infrared Contrast -- NEAR vs FAR
  
  
  c2all = where(emaf.c GE 0.2,nc2)
  c2far = where((emaf.c GE 0.2) AND (far OR tan),nc2f)
  
  print,'NC2: ',nc2,'  NC2F: ',nc2f
  
  c1all = where(emaf.c LT 0.1,nc1)
  c1far = where((emaf.c LT 0.1) AND (far OR tan),nc1f)
  
  print,'NC1: ',nc1,'  NC1F: ',nc1f
  
  bin = 0.02
  sbin = string(bin,format="(F0.2)")

  plothist,emaf[nind].c,bin=bin,color=ncolor,xr=[0,0.75],$
           xtit='Mid-Infrared Contrast',ytit='N per '+sbin+' contrast bin',$
           charsize=cs,thick=3,yr=[0,65],/yst,yminor=2,xmargin=[7,1],$
           ymargin=[4,1]
  
  ;; Far sources, filled histogram w/ median line
  plothist,emaf[find].c,bin=bin,color=fcolor,thick=3,/overplot,/fill,$
           fcolor=fcolor
  vline,median(emaf[find].c),linestyle=4,thick=3,color=fcolor

  ;; Near sources, outline histogram w/ median line
  plothist,emaf[nind].c,bin=bin,color=ncolor,thick=3,/overplot
  vline,median(emaf[nind].c),linestyle=3,thick=3,color=ncolor
  
  ;; Replot bottom x-axis to be seen atop the filled histogram
  cgAxis,xaxis=0,xthick=3,xr=[0,0.75],/xst,xtickformat='blank_axis'
  ;; axis,yaxis=0,ythick=3,yr=[0,70],/yst,ytickformat='blank_axis'
  
  print,'CONTRAST: ',median(emaf[nind].c),median(emaf[find].c)
  
  ;; Legend and panel label
  al_legend,/top,/right,box=0,linestyle=[0,0],thick=4,color=[ncolor,fcolor],$
            ['Near','Far + Tan'],charsize=cs*0.85,linsize=linsize
  al_legend,/top,/left,box=0,['(c)'],charsize=cs,/clear




  myps,/done
  
  












  ;; ;; Other junk




  ;; ;;====================================
  ;; ;; Panel 7: BGPS flux vs d_ML PLUS CONTRAST
  ;; myps,'./emaf_paper/plots/source_hists4.eps',xsize=10
  ;; multiplot,[2,1],xgap=0.035

  
  ;; far  = dml[*,0] GT tandist+0.5
  ;; near = dml[*,0] LT tandist-0.5
  ;; tan  = ~far AND ~near
  
  ;; print,'TOTAL: ',n_elements(far)
  ;; print,'NEAR:  ',fix(total(near))
  ;; print,'FAR:   ',fix(total(far))
  ;; print,'TAN:   ',fix(total(tan))
  
  ;; cgPlot,emaf.s_ap,dml[*,0],psym=16,symsize=0.5,ytit='d!dML!n [kpc]',/xlog,$
  ;;        xtit='BGPS 40" Flux Density [Jy]',charsize=1.0,/nodata,xr=logr,/xst
  
  
  ;; cgLoadct,40,/silent,/reverse
  ;; color = fix(round(emaf.c / 0.5 * 254.)<254)>0
  
  ;; FOR ij=0,n_elements(emaf)-1 DO $
  ;;     cgOplot,emaf[ij].s_ap,dml[ij,0],psym=16,symsize=0.5,color=color[ij]
  
  ;; cgColorbar,/vert,/right,charsize=0.6,divis=0,range=[0,0.5],$
  ;;            title='Contrast',position=[0.96,!Y.window[0],0.97,!Y.window[1]]

  ;; multiplot,/doyaxis
  
  ;; dx /= 10.
  ;; nx = n_elements(xarr)*10 - 1
  ;; xarr = findgen(nx)*dx + min(xarr)
  
  
  ;; s11 = 10.^(xarr)
  ;; d11 = findgen(149)+0.1
  ;; print,m4_stat(d11)
  
  ;; sarr = cmreplicate(s11,n_elements(d11))
  ;; darr = transpose(cmreplicate(d11,n_elements(s11)))
  
  ;; cube = readfits('./local/MW_model_ffore.fits',hd)
  ;; ff = cube[180,19,*]
  ;; farr = transpose(cmreplicate(ff[0:148],n_elements(s11)))
  ;; cont20 = (1.d - farr)*(1. - exp(-flux2tau_bgps(20.d)*sarr))

  ;; ff = cube[140,19,*]
  ;; farr = transpose(cmreplicate(ff[0:148],n_elements(s11)))
  ;; cont30 = (1.d - farr)*(1. - exp(-flux2tau_bgps(20.d)*sarr))

  ;; ff = cube[100,19,*]
  ;; farr = transpose(cmreplicate(ff[0:148],n_elements(s11)))
  ;; cont40 = (1.d - farr)*(1. - exp(-flux2tau_bgps(20.d)*sarr))


  
  ;; plotimage,cont30,range=[0.,0.5],xst=4,yst=4,min_dpi=100
  
  ;; cgContour,cont20,levels=findgen(5)/10.+0.1,/over,$
  ;;           c_color=[replicate('opposite',3),replicate('bisque',2)]

  ;; ;; cgContour,cont40,levels=findgen(5)/10.+0.1,/over,$
  ;; ;;           c_color=[replicate('opposite',3),replicate('bisque',2)]
  
  ;; cgAxis,charsize=1.0,xaxis=0,xr=logr,/xlog,/xst,$
  ;;        xtit='BGPS 40" Flux Density [Jy]',/save
  ;; cgAxis,xaxis=1,xtickformat='blank_axis'
  ;; cgAxis,yaxis=0,charsize=1.0,yr=[0,15],/save,$
  ;;        ytit='d!dML!n [kpc]'
  ;; cgAxis,yaxis=1,ytickformat='blank_axis'
  
  
  
  ;; myps,/done,/mp
  

  ;; ;;====================================
  ;; ;; Panel 9: BGPS flux vs d_ML PLUS CONTRAST PLUS MASS!
  ;; myps,'./emaf_paper/plots/source_hists5.eps',xsize=10
  ;; multiplot,[2,1],xgap=0.035

  
  ;; far  = dml[*,0] GT tandist+0.5
  ;; near = dml[*,0] LT tandist-0.5
  ;; tan  = ~far AND ~near
  
  ;; cgPlot,emaf.s_ap,dml[*,0],psym=16,symsize=0.5,ytit='d!dML!n [kpc]',/xlog,$
  ;;        xtit='BGPS 40" Flux Density [Jy]',charsize=1.0,/nodata,xr=logr,/xst
  
  
  ;; cgLoadct,40,/silent,/reverse
  ;; color = fix(round(emaf.c / 0.5 * 254.)<254)>0
  
  ;; FOR ij=0,n_elements(emaf)-1 DO $
  ;;     cgOplot,emaf[ij].s_ap,dml[ij,0],psym=16,symsize=0.5,color=color[ij]
  
  ;; cgColorbar,/vert,/right,charsize=0.6,divis=0,range=[0,0.5],$
  ;;            title='Contrast',position=[0.96,!Y.window[0],0.97,!Y.window[1]]

  ;; multiplot,/doyaxis
  
  ;; m = 13.1 * emaf.s_ap * dml[*,0]*dml[*,0]
  
  ;; cgPlot,charsize=1.0,m,dml[*,0],/xlog,ytit='d!dML!n [kpc]',/nodata,$
  ;;        xtit='M!dML!n [M!d'+cgSymbol('sun')+']',xtickformat='exponent10'
  
  ;; FOR ij=0,n_elements(emaf)-1 DO $
  ;;     cgOplot,m[ij],dml[ij,0],psym=16,symsize=0.5,color=color[ij]
  
  
  
  ;; myps,/done,/mp
  
END
